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Abstract 

Background: Competing risks methodology allows for an event-specific analysis of the single components of 
composite time-to-event endpoints. A key feature of competing risks is that there are as many hazards as there are 
competing risks. This is not always well accounted for in the applied literature. 

Methods: We advocate a simulation point of view for understanding competing risks. The hazards are envisaged 
as momentary event forces. They jointly determine the event time. Their relative magnitude determines the event 
type. 'Empirical simulations' using data from a recent study on cardiovascular events in diabetes patients illustrate 
subsequent interpretation. The method avoids concerns on identifiability and plausibility known from the latent 
failure time approach. 

Results: The 'empirical simulations' served as a proof of concept. Additionally manipulating baseline hazards and 
treatment effects illustrated both scenarios that require greater care for interpretation and how the simulation 
point of view aids the interpretation. The simulation algorithm applied to real data also provides for a general tool 
for study planning. 

Conclusions: There are as many hazards as there are competing risks. All of them should be analysed. This 
includes estimation of baseline hazards. Study planning must equally account for these aspects. 



Background 

The analysis of time -to -event data ('survival analysis') 
has evolved into a well established application of 
advanced statistical methodology in medicine. E.g., in 
the New England Journal of Medicine, survival methods 
have evolved from an occasionally used technique in the 
late 70s over moderate use in the late 80s into the lead- 
ing statistical procedure by 2005 [1]. The archetypical 
application analyses time until death, but combined end- 
points are also frequently considered. E.g., a recent lit- 
erature review in clinical oncology [2] found a multitude 
of combined endpoints including, e.g., progression-free 
survival, distant metastasis-free survival, locoregional 
relapse-free survival, etc. The medical problems at hand 
will, as these endpoints exemplarily suggest, usually be 
more complex than can be addressed by the analysis of 
time until one potentially combined event type. 
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Competing risks techniques allow for a more specific 
analysis in that they consider time until occurrence of 
the combined endpoint and endpoint type, e.g., progres- 
sion in contrast to death without prior progression. The 
relevance of competing risks in medical research is high- 
lighted by methodological papers in various medical 
fields. We mention [3-5] as recent examples. A classical 
statistics textbook account has been given in the first 
edition of [6] in 1980, and a definite mathematical treat- 
ment based on counting processes is included in [7]. 
Excellent tutorial papers in the statistical literature are 
[8,9]. 

However, despite an obvious practical relevance and a 
firmly established methodological background, compet- 
ing risks are not always well accounted for in published 
survival analyses in medical journals: E.g., another recent 
literature review [10] in clinical oncology found that 27 
out of 125 included randomised controlled trials consid- 
ered time to progression, but only 5 out of 125 studies 
accounted for such endpoint types being non-exclusive. 
A similar picture has been reported for studies on the 
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effect of implantable cardioverter defibrillator on subse- 
quent cardiac events [11]. 

The key to an adequate competing risks analysis is: 
There are as many hazards as there are competing risks. 
Unless all of these hazards, which are often called 
cause-specific hazards, have been analysed, the analysis 
will remain incomplete. In particular, only a complete 
analysis will allow for predicting event probabilities. 

The aim of this paper is to suggest an algorithmic or 
simulation point of view towards this key issue. The idea 
behind this point of view is that it gives a clear building 
plan of how the involved hazards generate competing 
risks data over the course of time. However, although the 
algorithm is mathematically well established [12], it is 
rarely used in practice [13]. We will use it as an opera- 
tional device for understanding competing risks. 

The remainder of the paper is organised as follows: The 
Methods Section introduces competing risks as arising 
from a multistate model, gives an algorithmic interpreta- 
tion of the involved hazards and provides a description of 
the data example from a recent study on cardiovascular 
events in diabetes patients [14]. In this Section, we also 
explain 'empirical simulation', where one simulates based 
on the empirical study probabilities, and give an overview 
on simulation scenarios that will be considered. The Sec- 
tion closes with a brief summary of the difference 
between the present paper and the common latent failure 
time approach to competing risks. Results are reported in 
the Results Section. Our simulations are 'empirical simu- 
lations', and they work as a proof of concept: Interpreting 
a competing risks analysis from a simulation perspective 
can be viewed as a thought experiment. The actual simu- 
lations show that the simulation perspective works. 
Finally, a discussion and a conclusion are offered, with a 
focus on consequences for practical competing risks ana- 
lyses, including graphical presentation and planning in 
the presence of competing events. 

Methods 

The competing risks multistate model 

Figure 1 depicts the multistate model of a competing 
risks process {X t ) t2 o with initial state 0 and two 
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Figure 1 Competing risks multistate model Competing risks 
multistate model with cause-specific hazards a 0j {f),j =1,2. 



competing event states 1 and 2. X t denotes the state 
that an individual is in at time t. The restriction to two 
competing events is for ease of presentation only. Initi- 
ally, every individual is in state 0 at time origin, i.e., X 0 
= 0. The individual stays in this state, i.e., X t = 0 until 
occurrence of any first event. Usually, there is one event 
of interest, modelled by transitions into state 1, and all 
other first event types may be subsumed into the com- 
peting event state 2. 

The competing risks process moves out of the initial 
state 0 at the event time T. At time T , the process 
enters one of the competing event states 1 or 2. Hence, 
the state X T that an individual is in at time T denotes 
the type of the first event, often called cause of failure. 
X T equals either 1 or 2. 

Key quantities in competing risks are the cause-speci- 
fic hazards (CSHs) a oj (t), j = 1, 2, 



aoj(t)dt = P(T e dt,X T = j, \T > i), j = 1, 2. 



(1) 



The dt notation in (1) is a short, but more intuitive 
form of the more common a 0 j(t) = lim Ats , 0 P(Ie [t, t + 
At), X T = / | T > t)/At. The interpretation of (1) is that 
the CSH at time t times the length dt of a very (infinite- 
simally) small time interval equals the conditional prob- 
ability of making an 0 — > / transition within that very 
small time interval [t, t + dt). The CSHs can be envi- 
saged as forces of transition, which work along the 
arrows in Figure 1. Analogous to cumulative distribution 
functions, basic statistical inference addresses cumula- 
tive quantities. We therefore also write A 0 j (t) for the 
cumulative CSHs, A oj (t) = A Q j(t) = f 0 a 0 j(u)du, j = 1,2. 

The CSHs are key, because, as seen below, the perti- 
nent probabilities are deterministic functions of the 
CSHs, and because both nonparametric estimation and 
Cox-type regression modelling build on them. The fun- 
damental nonparametric estimator is the Nelson- Aalen 
estimator of the cumulative CSHs 
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/ = 1, 2, where the summation is over all observed 
event times s < t. Variance estimation and, properly 
standardised, asymptotic normality of Aqj is discussed in 
detail in Chapter IV of [7], The most important regres- 
sion approach are Cox models for the CSHs. If, as is 
usual, we assume that covariates have different, i.e., 
cause-specific effects on the CSHs, the partial likelihood 
factorises such that a Cox model for a 0 i{t) m ^y be fitted 
by additionally censoring type 2 events, and vice versa 
for a Cox model for a 02 (t), see, e.g., Chapter VII of [7]. 
The R [15] packages mvna[16] and survival[17] 
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provide for convenient computation of the Nelson- Aalen 
estimator and the Cox model, respectively. 

The a 0 /s sum up to the all-cause hazard (Xo .(t)dt = P 
(fe df | T > t) with cumulative all-cause hazard A 0 .(t). 
The survival function of T is therefore a function of 
both a 0 /s, P(T > t) = exp(- fgU 0 .{u)du). As usual, P (T 
>t) is estimated using the Kaplan-Meier estimator, 
which is a deterministic function of the cause-specific 
Nelson-Aalen estimators. 

The so-called cumulative incidence functions (CIFs) 
are the expected proportion of individuals experiencing 
a certain competing event over the course of time, 

P(T < t,X T =))= f P(T > u-)a 0j (u)du, (2) 
Jo 

j = 1, 2. The interpretation of the right hand side of (2) is 
that the CIF at time t equals the integral over all 'probabil- 
ities' to make an 0 — > j transition precisely at time u. For 
this, an individual has to stay in state 0 until just before 
time u, which is reflected by the term P (T >u-). Condi- 
tional on being in state 0 just before time u, an 0 — > /' tran- 
sition at time u is reflected by the term a 0 j (u) du. 

One consequence is that the CIFs are involved functions 
of all CSHs via P (T >u-). The CIFs P{T<t,X T = 1) and 
P (T < t, X T = 2) add up to the all-cause distribution func- 
tion 1 - P (T >t). The Aalen-Johansen estimators of the 
CIFs can be obtained from (2) by substituting P {T >u-) 
with the Kaplan-Meier estimator and a 0 j (u) du with the 
increment of the cause-specific Nelson-Aalen estimator. A 
detailed discussion of the Aalen-Johansen estimator is in 
Chapter IV of [7]; variance estimation is assessed in [18]. 
In R, survival is typically used for estimating P (T >t). 
The Aalen-Johansen estimator may be computed using 
etm[19], and prediction of the CIFs based on Cox models 
for the CSHs is implemented in mstate[20]. 

It is crucial to any competing risks analysis that both 
cause-specific hazards completely determine the sto- 
chastic behaviour of the competing risks process [7, 
Chapter II. 6]. This is mirrored above by the CIFs and 
the all-cause distribution function being deterministic 
functions of all CSHs. However, these deterministic rela- 
tionships are involved. Therefore, we will pursue an 
algorithmic interpretation of the CSHs below. 

Algorithmic interpretation of the cause-specific hazards 

Thinking of the CSHs as momentary forces of transition 
which move along the arrows in Figure 1 suggests that 
competing risks data are generated over the course of 
time as follows: 

1. The event time T is generated with distribution 
function 1 - P (T >t), i.e., with hazard a 01 (t) + a 02 (t) 
= a 0 .(t). 



2. At time T , event type /' occurs with probability 
ay (T)/a 0 .(T),j = 1, 2. 

Using this algorithm for simulation studies in compet- 
ing risks has been discussed in [13]. It is important to 
note, however, that the algorithm goes beyond the com- 
putational question of how to implement simulations. 
Rather, the algorithm reflects the probabilistic question 
of how to build a probability measure based on the 
CSHs. This aspect is discussed in detail by [12] in the 
more general context of multistate models which are 
realized as a nested series of competing risks 
experiments. 

The algorithmic perspective of this paper then implies 
that the task of statistical inference is to detect the 
ingredients of the above algorithm. We illustrate this 
approach in the data example below. To this end, we 
note that the analysis of a combined endpoint is 
restricted to step 1 of the above algorithm. Here, the 
effect of a treatment, say, on both CSHs determines 
whether the occurrence of an event (of any type) is 
delayed or accelerated. In step 2, the type of an event 
again depends on the treatment effect on both CSHs. 
We illustrate below that interpretation is straightforward 
if the treatment effects on the CSHs work in opposite 
directions or if one CSHs remains unaffected. However, 
interpretation will become more challenging if there are 
unidirectional effects on both CSHs. We will find that, 
in general, it is also mandatory to consider cause-speci- 
fic baseline hazards in the interpretation. 

The 4D study 

The background of the 4D study was that statins are 
known to be protective with respect to cardiovascular 
events for persons with type 2 diabetes mellitus without 
kidney disease, but that a potential benefit of statins in 
patients receiving hemodialysis had until then not been 
assessed. Patients undergoing hemodialysis are at high 
risk for cardiovascular events. The 4D study was a pro- 
spective randomised controlled trial evaluating the effect 
of lipid lowering with atorvastatin in 1255 diabetic 
patients receiving hemodialysis. Patients with type 2 dia- 
betes mellitus, age 18-80 years, and on hemodialysis for 
less than 2 years were enrolled between March 1998 
and October 2002. Patients were randomly assigned to 
double-blinded treatment with either atorvastatin (619 
patients) or placebo (636 patients) and were followed 
until death, loss to follow-up, or end of the study in 
March 2004. 

The 4D study was planned [21] and analysed [14] for 
an event of interest in the presence of competing risks. 
The event of interest was defined as a composite of 
death from cardiac causes, stroke and non-fatal myocar- 
dial infarction, whichever occurred first. The other 
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competing event was death from other causes. Wanner 
et al. reported a CSH ratio of 0.92 (95%-confidence 
interval [0.77, 1.10]) for the event of interest. There was 
essentially no effect on the competing CSH. 

The simulation study below will use the data of the 
placebo group. In this group, 243 (38.2% of 636) events 
of interest and 129 (20.3%) competing events were 
observed. There were 264 (41.5%) censored patients. 
The Nelson-Aalen estimators Aoi(t) f° r the event of 
interest and Ao2(t) 101 the competing event are dis- 
played in Figure 2. 

'Empirical simulation' 

Simulations were based on the empirical probabilities 
defined by (Ani(t), Ao2(t)} ^° avo ^ the issue of model 
misspecification, which is outside the scope of the pre- 
sent investigations, 'ideal' Cox models were used for 
generating data in the treatment group as explained 
below. Event times and event types were generated fol- 
lowing the two-step algorithm described earlier. 

To be specific, event times in the placebo group were 
drawn from the distribution defined by the Kaplan- 
Meier estimator derived from (Aoi(t), Ao2(t))- Here, a 



practical complication arose in that the Kaplan-Meier 
estimator only spent about 76% of the probability mass. 
I.e., the Kaplan-Meier curve did not drop down to 0 
because of right-censoring, which is a common phe- 
nomenon in clinical trials. This was handled by putting 
point mass beyond the largest observed time; corre- 
sponding realisations were always censored. In other 
words, on average about 100% - 76% = 24% of the simu- 
lated event times equalled some time that was larger 
than the largest observed time in the original data set. 
The corresponding observations were always censored 
with censoring times generated as explained below. 
Event types were generated by the binomial experiment 
of the two-step algorithm, substituting the CSHs by the 
increments of the cumulative CSHs. The algorithm was 
analogously applied, when transformed cumulative 
CSHs were used for the placebo group. 

Data in the treatment group were generated based on 
'ideal' Cox models. That is, CSH ratios exp(/3i) for the 
event of interest and exp(/J 2 ) for the competing event 
were specified. If the original empirical baseline hazards 
were used, data were drawn based on the probabilities 

defined by (exp(fr)Aoi(t), exp^A^f)} If 
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Figure 2 Nelson-Aalen estimators of the cumulative CSHs. Nelson-Aalen estimators (solid lines) of the cumulative CSHs with log-transformed 
95% confidence intervals in the placebo group. 
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transformations of (Aoi(t), Arj2(f)) were used for gener- 
ating placebo data, the CSH ratios acted on the trans- 
formed cause-specific baseline hazards. This reflects the 
situation that the true underlying baseline CSHs are a 0 i 
(t) and a 02 {f) (or transformations thereof) in the control 
group, while the CSHs of the treatment group are exp 
(Pi)a 0 i{t) and exp(P 2 )a 02 (t) for the case of untrans- 
formed baseline hazards. Random censoring times were 
generated for all individuals based on the Kaplan-Meier 
estimator of the censoring survival function in the pla- 
cebo group. Note that this estimator did spend 100% of 
the probability mass. 

Overview of simulation scenarios 

The scenarios investigated in the next Section differed 
both with respect to the choice of Pi and P 2 and in 
terms of the cumulative baseline hazards, i.e., the cumu- 
lative CSHs in the placebo group. The choice of the P's 
broadly falls into three categories. 

One category is characterised by p 2 = 0, i.e., there is 
no effect on the competing CSH. A prime example are 
implantable cardioverter defibrillators [11], which dis- 
play a beneficial effect on the CSH of sudden cardiac 
death but no effect on the CSH for death from other 
causes. The assumption of p 2 = 0 straightforwardly 
implies the direction of the treatment effect on the CIF: 
If Pi < 0, the CIF of interest in the treatment group is 
always less than the one in placebo group. The relation- 
ship is reversed, if Pi > 0. This is intuitively understood 
thinking of the CSHs as momentary forces of transition, 
and it is reflected in step 2 of the simulation algorithm. 
E.g., if Pi < 0 and P 2 = 0, the binomial event type 1 
probabilities are reduced for the treatment group. 

A second category is characterised by opposite treat- 
ment effects on the CSHs. This category straightfor- 
wardly implies the direction of the treatment effect on 
the CIF, too: The constellation Pi < 0 and P 2 > 0 implies 
a smaller CIF of interest in the treatment group but also 
a larger competing CIF. These relations are reversed for 
Pi > 0 and P 2 < 0. This is again intuitively implied by 
thinking of the CSHs as momentary forces of transition, 
and it is also reflected in step 2 of the simulation algo- 
rithm. E.g., if Pi < 0 and P 2 > 0, the binomial event type 
1 probabilities are reduced for the treatment group. 

Finally, unidirectional treatment effects on the CSHs 
constitute the third category. Interestingly, the interpre- 
tation of unidirectional effects is straightforward when 
both competing events are fatal in the sense that a treat- 
ment with Pi < 0 and p 2 < 0, say, is beneficial. But uni- 
directional effects also present the most challenging 
scenario in terms of understanding the resulting course 
of the CIFs. The interpretation for step 1 of the simula- 
tion algorithm is straightforward. If, e.g., Pi < 0 and /J 2 
< 0, events of any type will happen later. The 



interpretational challenge becomes apparent in the sec- 
ond step of the algorithm: If, e.g., P 2 <Pi < 0, the relative 
magnitude of the CSHs changes such that the binomial 
event of interest probabilities are increased. The inter- 
pretational difficulty is the increase of this probability, 
although pi is negative. This constellation may result in 
an (eventually) increased CIF of interest. 

All three categories are encountered in practice. 
Examples from clinical trials in hospital epidemiology 
are given in [22]. 

Difference to the latent failure time model 

Some readers may be more familiar with competing risks 
as arising from risk-specific latent times, say, I™ and I< 2) . 
The connection to our multistate framework is T = min(X 
r< 2) ) with event type X T = 1, if T* 1 ' <J (2 \ and X T = 2 
otherwise. The latent failure time model imposes an addi- 
tional structure, which has been heavily criticised mainly 
for three reasons. The dependence structure of and T 
^ is, in general, not identifiable [23]. Since the latent times 
are unobservable, there is something hypothetical about 
them, which questions their plausibility [24] . Perhaps most 
importantly, it has been disputed whether the latent failure 
time point of view constitutes a fruitful approach to answer 
questions of the original subject matter [25]. 

Despite of this critique, latent times are the predomi- 
nant approach for simulating competing risks data [13]. 
Assuming, for tractability, T- and 1^ to be indepen- 
dent with hazards equal to the cause-specific hazards 
aoi(i) and (Xo 2 (t), respectively, is computationally correct 
in that simulation based on this model yields the right 
data structure. 

However, nothing is gained from assuming the addi- 
tional latent structure, either. As a consequence, we will 
emphasise simulation and interpretation along the lines 
of the Algorithmic interpretation of the cause-specific 
hazards outlined earlier. This avoids the concerns on 
identifiability, plausibility and usefulness. 

E.g., in the 4D study, the typical interpretation of the 
latent times would be that is the time until death from 
cardiac causes, stroke or non-fatal myocardial infarction, 
while is the time until death from other causes. Such 
an interpretation has given rise to debating whether, say, a 
patient may conceptually still die from other causes after 
having died because of a cardiac event. In contrast to this, 
our approach only assumes that a patient is conceptually 
at risk of experiencing any of these events, provided that 
none of these events has happened so far. 

Results 

General 

We studied ten different scenarios, which are tabulated 
in Table 1 and cover all effect categories discussed ear- 
lier. Scenarios 1-5 have Pi similar to the actual study 
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Table 1 Scenarios in the simulation study 
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Scenarios in the simulation study. Scenarios 1-5 have /Ji similar to the actual 
study result, scenarios 6-10 have similar to the planning figures. 



result [14], and scenarios 6-10 have Pi similar to the 
planning figures [21]. Table 1 also displays the use of 
different baseline hazards. For each of the scenarios, 
1000 simulation runs were considered with 500 indivi- 
duals in the placebo group and 500 individuals in the 
treatment group. 

Table 2 presents, for each scenario and based on Cox 
analyses of the CSHs, the mean log CSH ratios p and 
j$ , empirical 95% confidence intervals, the coverage 
probabilities of the Wald-type confidence intervals for 
the estimated regression coefficients, and the empirical 
power. 

For Scenarios 1-5, we also plotted the true CIFs (solid 
black lines), the average of the Aalen-Johansen estima- 
tors (dashed black lines) and 300 randomly selected 
Aalen-Johansen estimators (grey lines). The reason to 
only plot a random subsample was to maintain a grey 
shading in the plots. 

We collect some general observations. In the Figures, 
the true CIFs are visually hardly distinguishable from 



the average of the Aalen-Johansen estimators from each 
simulation study. This entails that the algorithmic point 
of view is appropriate: Simulation along this line yields 
data that are consistent with the original quantities. 
Similar statements hold for the average of the estimated 
regression coefficients and the coverage probability of 
their confidence intervals. The Figures also show that 
both regression coefficients and baseline hazards matter. 
E.g., keeping both regression coefficients, but changing 
baseline hazards alters the CIFs. Similarly, keeping both 
the CSH ratio of interest and the baseline hazards, but 
changing the competing CSH ratio alters the CIFs. 

We also note that recovering the true CIFs implies 
that we would have also recovered the original cumula- 
tive hazards of Figure 2. This is so, because knowledge 
of all CIFs allows to derive all CSHs and vice versa. 

Scenarios 1 and 2: no effect on the competing CSH 

Scenarios 1 and 2 are chosen with regression coeffi- 
cients similar to the 4D study. Scenario 1 used the origi- 
nal cumulative baseline hazards (Arji(f), Arj2(t)) an d 
scenario 2 used the transformed tuple 
((Aoi(O) 1 ' 4 ' Ao2(t))' This transformation amplifies the 
hazard for the event of interest, because Ani(t) < 1 
except for the right tail of the time interval displayed in 
Figure 2. 

Figures 3 and 4 illustrate that the CIF of interest is 
lower in the treatment group, as implied by Pi < 0 
under the side condition p 2 = 0. However, the difference 
is small, because the effect mirrored by Pi = -0.1 is 
rather moderate. The effect is somewhat amplified for 
the CIFs of interest when using (Aoi(t)) 1 '' 4 a s a baseline 
CSH, because the transformation amplifies this hazard. 
The effect of a 'more important' baseline CSH for the 
event of interest is also reflected in the smaller empirical 
95% confidence interval for the estimation of Pi and an 



Table 2 Simulation results 
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-0.31; 0.11 


94.9 


15.1 
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0.03; 0.56 


94.5 


64.5 
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-0.1 


-0.3; 0.09 


94.8 


16.4 


0.3 


-0.09; 0.73 


96.2 


27 
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-0.1 


-0.38; 0.21 


96.1 


11.1 


-0.31 
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93.7 


91.4 
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94.1 


78.8 
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94.5 


95.1 
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-0.29 


-0.51; -0.09 


95.7 


75.9 


0.29 


0.04; 0.56 
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59.2 
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-0.3 


-0.48; -0.09 


95.5 


81.4 


0.31 


-0.11; 0.78 


94.5 


26.9 


10 


-0.3 


-0.61; 0.03 


94.4 


46.9 


-0.32 


-0.49; -0.14 
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92.1 



Simulation results. For each scenario and for each competing risk are presented the averaged log hazard-ratios, empirical confidence intervals (CI), coverage 
probabilities and empirical power. 
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increased empirical power. This is so, because increasing 
the baseline CSH of interest while not changing the 
competing baseline CSH will lead to more events of 
interest. In both scenarios, we find a slight increase of 
the competing CIF in the treatment group. Comparing 
both scenarios, one also finds that the overall magnitude 
of the competing CIF is reduced by amplifying the base- 
line CSH of interest. Finally, we note that the empirical 
power for the competing event approximately keeps the 
nominal level of 0.05. 

Scenarios 3 and 4: opposite treatment effects on the 
CSHs 

Scenarios 3 and 4 are chosen with /?i similar to the 4D 
study, but /3 2 > 0 having an opposite effect. We deliberately 
chose \p 2 \ > \Pi\ and, in scenario 4, transformed cumula- 
tive baseline hazards (Ani(t), (Ao2(£)) 2 ) to illustrate that 
the magnitude of a regression coefficient must be seen in 
connection with the magnitude of the corresponding base- 
line hazard. The transformation x >-* x reduces the magni- 
tude of the competing hazard, because Arj2(t) < 1- Figures 
5 and 6 confirm a decreasing treatment effect on the CIF 
of interest and an increasing effect on the competing CIF. 
These effects are, however, diminished when using 
(Arj2(0) 2 as a baseline CSH, because the transformation 
reduces the magnitude of this hazard. This is, e.g., 



reflected in the empirical 95% confidence interval for the 
estimation of /? 2 , which excludes 0 only in the case of 
untransformed baseline CSHs. In analogy to this, the 
empirical power is considerably reduced for the competing 
event in the transformed case. 

Scenario 5: unidirectional treatment effects on the CSHs 

We chose fii = -0.1 as before and f} 2 = -0.3. Figure 7 
displays the true CIFs with transformed cumulative 
baseline hazards ((Aoi(t)) 2 , (Ao2(t))^ 4 )- Because the 
transformation amplifies the competing CSH, but has an 
opposite effect on the CSH of interest, we see a slowed 
down increase of the CIF of interest in the control 
group as compared to the first two scenarios. There is a 
visible treatment effect on the competing CIF, which is 
reduced as compared to the control group. In contrast, 
the CIF of interest evolves at a comparable and low 
magnitude in both groups before eventually displaying 
larger probabilities for the treated. In fact, the CIFs of 
interest cross, i.e., the curve for the treated first runs 
below the one for the control group and then crosses. 
However, the early difference is hardly visible due to the 
overall low magnitude of the CIF of interest. 

The eventual difference between the CIFs of interest 
appears to contradict Pi < 0, but can well be understood 
from a simulation perspective as outlined in the 



Interest - Treated 




Competing - Treated 




Interest - Untreated 




Competing - Untreated 




Figure 5 Scenario 3 of Table 1, Solid black lines are the true CIFs, solid grey lines are 300 randomly selected Aalen-Johansen estimators. The 
average of the Aalen-Johansen estimators is drawn as dashed black lines. Dotted lines in the left plots are the corresponding untreated CIFs. 
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Overview of simulation scenarios. This is illustrated in 
Figure 8. We discuss the left part of the Figure first. It 
corresponds to the first step of the simulation algorithm 
and illustrates that events of any type happen later in 
the treatment group, because both /3 1 and /3 2 are 
negative. 

The plot shows the point masses of the empirical 
event time distribution for both treatment groups. The 
Figure has been restricted to [0.07, 5.4] in order to 
avoid the steep increases of the survival function 
towards the corners of the observed time interval. Black 
bars indicate larger point mass as compared to their 
corresponding bars in the other group. The correspond- 
ing bars indicating smaller point mass are grey. We also 
note that black bars do not really superimpose grey bars 
and vice versa in the Figure, although this impression is 
occasionally conveyed by the high density of bars. Figure 
8 (left) clearly shows that probability mass is moved 
towards later event times in the treatment group. This 
is also reflected by the CIFs, which start to increase 



later for the treatment group, although this is hardly 
visible for the CIF of interest. 

However, the left plot does not illustrate why the CIFs 
of interest cross. This is explained in the right plot. It 
corresponds to the second step of the simulation algo- 
rithm and illustrates that treatment shifts probability 
mass towards type 1 events in general and, more specifi- 
cally, to later type 1 events as a consequence of f} 2 <fii < 
0, see the Overview of simulation scenarios. 

To this end, note that in an actual data situation the 
empirical counterpart of step 2 of the simulation algo- 
rithm will often equal either 0 or 1, if there are no type 
1 and type 2 events being observed at the same time. In 
our data example, both type 1 and type 2 events hap- 
pened at only 18 out of 322 time points considered in 
Figure 8. Figure 8 (right) profits from the fact that the 
binomial probabilities are often degenerated and only 
shows those times with an observed type 1 event. If 
such a time is drawn, the event type is almost always 
determined to be of type 1. 



All even! times in [0.07, 5.4] shown 



Only times in [0.07, 5.4] with type 1 events shown 
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Figure 8 Scenario 5, point masses of the empirical event time distribution. Point masses in [0 07, 5.4] of the empirical event time 
distribution. The left plot shows all event times. The right plot is restricted to times with an event of interest. Bottom of the plots are for the 
control group with ordinate indicated on the left, top of the plots for the treatment group with ordinate indicated on the right. Black bars are 
larger than their corresponding bars in the other group, which then are grey. 
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The right plot illustrates two things: Firstly, black 
colour dominates the upper part of the plot, indicating 
that the event of interest is more likely to happen in 
the treatment group. Secondly, the colouring moves 
from grey to black for the treatment group and from 
black to grey in the control group. The interpretation 
is that, initially, event type 1 times are drawn with 
higher probability in the control group. The picture is 
being reversed as time proceeds, which leads to cross- 
ing CIFs, and eventually the cumulative proportion of 
type 1 events is larger in the treatment group. Note, 
however, that there is overall low probability mass on 
early type 1 event times, which implies that the CIFs 
of interest initially are hardly distinguishable and 
rather small. 

An analogous plot of Figure 8 (right) for type 2 events 
shows that probability mass is almost uniformly reduced 
for type 2 event times by the treatment effect (figure 
not shown). 

Scenarios 6-10 

Scenarios 6-10 repeated the previous investigations with 
a more pronounced treatment effect of f5\ = -0.3. Results 
are reported in Table 2. The most striking difference to 
the results from scenarios 1-5 is an increased empirical 
power for events of type 1. The increased empirical 
power, however, does not only depend on fi 1 = -0.3, but 
on all aspects discussed above. E.g., both scenarios 6 
and 7 have (/J 1( fi 2 ) = (-0.3, 0), but the cause-specific 
baseline hazard for type 1 events is amplified in scenario 
7. This leads to scenario 7 having better empirical 
power than scenario 6. In contrast, power is substan- 
tially decreased for the situation studied in scenario 10. 

Discussion 

This paper envisaged the CSHs as momentary forces of 
transition, which suggests an algorithmic perspective 
towards competing risks. 'Empirical simulations' worked 
as a proof of concept. The algorithmic perspective was 
used on the interplay between CSHs and CIFs. 

The involved relationship between CSHs and CIFs has 
inspired a boost in methodological research on testing 
and direct modelling of the CIFs, e.g., [26-29]. Similar to 
our paper, a number of recent references have used 
simulation of competing risks data to investigate these 
methods. In particular, Gray's [26] test [5,30,31] and the 
Fine-Gray [27] model [22,32] have attracted attention. 
Both these exemplary references and the present paper 
found a subtle interplay between different CSH constel- 
lations and subsequent impact on the CIFs. 

The difference to the present paper is that a typical 
simulation study will put the simulation algorithm aside 
as only a computational tool, once the data have been 
generated. Thus, one will typically specify the CSHs, use 



some simulation algorithm [13] for data generation and 
analyse the data with the methodology at hand. Then, 
the CSH specifications and the results of the data ana- 
lyses will be compared. In contrast to this, we have 
advocated to use the simulation algorithm itself as an 
operational tool for interpretation. In this context, it is 
worthwhile to note that the algorithm of our paper does 
not experience a number of problems which come with 
the common latent failure time model. E.g., the problem 
of dependence of the latent times has motivated to 
include different dependence structures in some simula- 
tion studies. There is no such problem in our set-up, 
which therefore facilitates interpretation. 

We discuss practical consequences next. To begin, it 
is interesting to revisit the results of the 4D study in the 
light of the simulation algorithm. As stated earlier, the 
original study finding was a CSH ratio of 0.92 for the 
event of interest and essentially no effect on the com- 
peting CSH. The competing CSH ratio was, of course, 
not exactly equal to 1.00, but it displayed a slight reduc- 
tion. Because the treatment effect on the CSH of inter- 
est was moderate, and because the Nelson-Aalen 
estimators of the cumulative CSHs were not exactly pro- 
portional between treatment groups, the Aalen-Johansen 
estimators of the CIFs of interest display a somewhat 
subtle relationship in the original report. (Figure three 
in [14], not reproduced here.) E.g., the CIFs cross before 
displaying a moderate benefit for the treatment group. 
However, the difference between the CIFs is slight 
before crossing and must therefore not be overinter- 
preted. As a consequence, we believe that interpretation 
of the competing risks situation at hand is well guided 
by the idealised situation of the simulation scenario 1. 

Next, we reiterate that it is crucial that all CSHs are 
analysed. In the Backgroung Section, we noted that this 
is often not the case in clinical research. We also illu- 
strated that a comprehensive analysis should not be 
restricted to hazard ratios only, but that ideally the 
cause-specific baseline hazards will be considered, too. 
The bottom line is that the interpretation of the CSH 
ratio of interest depends both on the baseline CSH of 
interest, the competing CSH ratio and the competing 
baseline CSH. In particular, a missing analysis of the 
competing CSH may have seriously misleading 
consequences. 

An important issue in this context is that of graphi- 
cally presenting results. A popular and adequate choice 
are plots of the CIFs, which should, in particular, be 
plotted for all event types, if all competing risks are 
harmful. However, it was also illustrated that the con- 
nection between these plots and CSH analyses is not 
straightforward, such that further graphical tools would 
be helpful. The most obvious choice is to also show the 
estimated cumulative CSHs as in Figure 2. This should 
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be done much more often. The reason is that the CSHs 
regulate the stochastic behaviour of the competing risks 
process as explained in the Methods Section and as illu- 
strated in the Results Section. 

In addition, plots such as Figure 8 which highlight the 
simulation perspective can be useful. The interested 
reader is also referred to 'vertical modelling' [33] and 
multistate incidence rate graphics [34]. We also note that 
'vertical modelling' aims at modelling the binomial prob- 
abilities in step 2 of the simulation algorithm as a smooth 
curve. This approach is appropriate when taking the 
simulation algorithm as a starting point to model com- 
peting risks data. We reiterate that our aim has been dif- 
ferent in that we took the simulation perspective as an 
operational tool to interpret the standard CSH analyses. 

We return to the key fact that all CSHs should be 
analysed in the presence of competing risks, and that 
this is often not accounted for in clinical research. 
These issues raise the question of planning competing 
risks studies, see [35] for a recent review. If the aim is 
to compare CSHs, a typical assumption made during the 
planning phase of the study is that of constant or piece- 
wise constant hazards, e.g., [21,36,37]. In addition, it is 
often assumed that the treatment does not affect the 
competing CSH, see [35]. In practice, it may be difficult 
to find an adequate closed form for time-dependent 
CSHs. Sample size calculations may become quite for- 
midable, if the planned analysis is more complex than 
testing the CSH of interest only. Whatever the planned 
statistical experiment is, we note that our empirical 
simulation approach provides for a general tool to study 
empirical power and, hence, to decide on sample size, if 
data of a control group - or of patients similar to the 
anticipated control group - are available. This is also 
illustrated in Table 2. Interestingly, Figure 2 suggests 
that assuming constant CSHs might be a reasonable 
assumption for the present control group, but it should 
be pointed out that the simulation approach does not 
rely on such an assumption. It could be applied without 
further ado, if the CSHs show a pronounced time- 
dependency. In particular, one would not need to spe- 
cify a closed form for the time-dependent CSHs. In clos- 
ing, we mention that, while we have focused on the Cox 
model as the major tool to analyse CSHs, other models 
such as Aalen's additive model may be used for CSHs, 
too; see [38] for a recent textbook treatment. The simu- 
lation perspective of this paper may then applied to 
results from other CSH models, too. 

Conclusions 

This paper suggests an algorithmic or simulation point 
of view for the interpretation of competing risks ana- 
lyses. This point of view follows the construction of 
competing risks data based on the CSHs, envisaging the 



hazards as momentary forces of transition. Concerns on 
identifiability and plausibility that are common in the 
latent failure time context do not arise. 

Simulation studies based on the empirical probability 
measure of a real data analysis served as a proof of con- 
cept. The simulation point of view was found to be ade- 
quate in that it recovered the original empirical law. 
Manipulating baseline hazards and treatment effects 
highlighted different aspects of a competing risks 
analysis. 

All CSHs should be analysed, including the cause-spe- 
cific baseline hazards. 'Empirical simulations' also pro- 
vide a flexible tool for study planning in the presence of 
competing risks. 
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